!
! Copyright (C) 2000-2013 A. Marini and the YAMBO team 
!              https://code.google.com/p/rocinante.org
! 
! This file is distributed under the terms of the GNU 
! General Public License. You can redistribute it and/or 
! modify it under the terms of the GNU General Public 
! License as published by the Free Software Foundation; 
! either version 2, or (at your option) any later version.
!
! This program is distributed in the hope that it will 
! be useful, but WITHOUT ANY WARRANTY; without even the 
! implied warranty of MERCHANTABILITY or FITNESS FOR A 
! PARTICULAR PURPOSE.  See the GNU General Public License 
! for more details.
!
! You should have received a copy of the GNU General Public 
! License along with this program; if not, write to the Free 
! Software Foundation, Inc., 59 Temple Place - Suite 330,Boston, 
! MA 02111-1307, USA or visit http://www.gnu.org/copyleft/gpl.txt.
!
integer function ioX(X,Xw,ID)
 !
 use pars,          ONLY:SP,schlen
 use units,         ONLY:HA2EV
 use memory_m,      ONLY:mem_est
 use R_lattice,     ONLY:nqibz,g_vec
 use matrix_operate,ONLY:mat_c2r,mat_r2c,mat_r2c_ordered
 use X_m,           ONLY:X_t,X_mat,X_RLcomp_ig,X_description
 use frequency,     ONLY:w_samp
 use IO_m,          ONLY:io_connect,io_disconnect,io_sec,io_com,&
&                        io_elemental,io_status,io_bulk,io_restart_point,&
&                        read_is_on,write_is_on,io_header,&
&                        RD_CL,RD_CL_IF_END,io_action,&
&                        ver_is_gt_or_eq,io_serial_number,serial_number,&
&                        NONE,IO_INCOMPATIBLE_VAR
 use stderr,        ONLY:intc
 use fragments,     ONLY:io_fragment,Fragments_Restart
 implicit none
 type(X_t)   :: X
 type(w_samp):: Xw
 integer     :: ID
 !
 !Work Space
 !
 integer             ::sec_size,i1,iq,i_err,io_com_save
 logical             ::different_db_RL_order
 character (schlen)  ::ch,db_desc
 character (2)       ::local_XC_KIND
 real(SP),allocatable::X_disk(:,:,:),W_disk(:,:),RL_vecs_disk(:,:)
 !
 ! I assume that he RL ordering used to build X is not different
 ! from the one in DB1
 !
 different_db_RL_order=.false.
 !
 if (X%whoami==1) db_desc='Xx'
 if (X%whoami==2) db_desc='em1s'
 if (X%whoami==3) db_desc='em1d'
 if (X%whoami==4) db_desc='pp'
 !
 ioX=io_connect(desc=trim(db_desc),type=2,ID=ID) 
 if (ioX/=0) goto 1
 !
 ! Check the restart point
 !
 call Fragments_Restart(ID)
 !
 if (any((/io_sec(ID,:)==1/))) then
   !
   ! Header
   !
   if (X%whoami==1) local_XC_KIND="Xx"
   if (X%whoami==2) local_XC_KIND="Xs"
   if (X%whoami==3) local_XC_KIND="Xd"
   if (X%whoami==4) local_XC_KIND="Xp"
   ioX=io_header(ID,QPTS=.true.,R_LATT=.true.,WF=.true.,T_EL=.true.,XC_KIND=local_XC_KIND,CUTOFF=.true.,GAUGE=.TRUE.)
   if (ioX/=0) goto 1
   !
   ! Global Descriptions
   !=====================
   !
   call variables_X(ID,X_description,X)
   !
   ! When the SN is different and I am importing this DB 
   ! from a different calculation I check the RL vectors
   !
   if (ver_is_gt_or_eq(ID,(/3,0,6/))) then
     !
     allocate(RL_vecs_disk(X%ng_db,3))
     !
     if (write_is_on(ID)) RL_vecs_disk(:,:)=g_vec(:X%ng_db,:)
     !
     call io_bulk(ID,VAR='X_RL_vecs',VAR_SZ=shape(RL_vecs_disk))
     call io_bulk(ID,R2=RL_vecs_disk)
     !
     if (io_serial_number(ID)/=serial_number.and.read_is_on(ID).and.&
&        allocated(g_vec)) then
       if (.not.allocated(X_RLcomp_ig)) allocate(X_RLcomp_ig(X%ng))
       i_err = X_RLcomp_table() 
       if ( i_err < 0 ) io_status(ID) = IO_INCOMPATIBLE_VAR
       if ( i_err > 0 ) different_db_RL_order=.true.
     endif
     deallocate(RL_vecs_disk)
     !
   endif
   !
   ! Close & Menu lower line
   !
   if (.not.any((/io_sec(ID,:)==2/))) then
     call io_elemental(ID,VAR="",VAR_SZ=0,MENU=1)
     ioX=io_status(ID)
     if (ioX/=0) goto 1
   else
     call io_elemental(ID,VAR="",VAR_SZ=0,MENU=0)
   endif
 endif
 !
 do i1=1,size(io_sec(ID,:))
   if (io_sec(ID,i1)==0) cycle
   iq=io_sec(ID,i1)/2
 enddo
 if (iq==0) goto 1
 !
 ! Manage RD_CL_IF_END
 !
 if (io_action(ID)==RD_CL_IF_END.and.iq==nqibz) io_action(ID)=RD_CL
 !
 ! Fragmentation
 !
 call io_fragment(ID,i_fragment=iq)
 !
 !  type(X_t)    ::X(4)  x s d p
 !
 if (any((/io_sec(ID,:)==2*iq/))) then
   sec_size=6
   !
   if (X%whoami==4) sec_size=1
   write (ch,'(a,i3.3)') "FREQ_PARS_sec_iq",iq
   if (ver_is_gt_or_eq(ID,revision=308)) ch="FREQ_PARS_sec_iq"//trim(intc(iq))
   !
   io_com_save=io_com(ID)
   if (io_sec(ID,1)/=1) io_com(ID)=NONE
   !
   if (io_sec(ID,1)==1) call io_elemental(ID,VAR=trim(ch),VAR_SZ=sec_size,MENU=0)
   if (io_sec(ID,1)/=1) call io_elemental(ID,VAR=trim(ch),VAR_SZ=sec_size,MENU=1)
   !
   if (X%whoami==4) then
     call io_elemental(ID,UNIT=HA2EV,&
&         VAR=" PPA Im  energy     [ev]:",R0=X%ppaE,CHECK=.true.,OP=(/"=="/))
   endif
   if (X%whoami/=4) then
     call io_elemental(ID,&
&         VAR=" :: Current Q-pt index     :",I0=iq)
     call io_elemental(ID,UNIT=HA2EV,&
&         VAR=" :: X energy range     [ev]:",R1=Xw%er,CHECK=.true.,OP=(/"==","=="/))
     call io_elemental(ID,UNIT=HA2EV,&
&         VAR=" :: X damping range    [ev]:",R1=Xw%dr,CHECK=.true.,OP=(/"==","=="/))
     call io_elemental(ID,&
&           VAR=" :: Number of frequencies  :",I0=Xw%n(1),CHECK=.true.,OP=(/"=="/))
   endif
   !
   call io_elemental(ID,VAR="",VAR_SZ=0,MENU=1)
   ioX=io_status(ID)
   !
   io_com(ID)=io_com_save
   !
   if (ioX/=0) goto 1
   !
 endif
 !
 ! Sections 1 & 2 are used in VERIFY mode. 
 ! Now that the menu is closed I can return with ioX/=0 if
 ! there is a restart point
 !
 ioX=io_restart_point(ID)
 if (ioX/=0) goto 1
 !
 if (any((/io_sec(ID,:)==2*iq+1/))) then
   allocate(X_disk(X%ng_db,X%ng_db,2))
   !
   write (ch,'(a,i3.3)') "FREQ_sec_iq",iq
   if (ver_is_gt_or_eq(ID,revision=308)) ch="FREQ_sec_iq"//trim(intc(iq))
   call io_bulk(ID,VAR=trim(ch),VAR_SZ=(/Xw%n(1),2/))
   allocate(W_disk(Xw%n(1),2))
   !
   if (write_is_on(ID)) then
     W_disk(:,1)=real(Xw%p(:))
     W_disk(:,2)=aimag(Xw%p(:))
   endif
   !
   call io_bulk(ID,R2=W_disk)
   !
   if (read_is_on(ID)) then
     if (associated(Xw%p)) then
       deallocate(Xw%p)
       call mem_est("W-p")
     endif
     allocate(Xw%p(Xw%n(1)))
     call mem_est("W-p",(/Xw%n(1)/))
     Xw%p(:)=cmplx(W_disk(:,1),W_disk(:,2))
   endif
   !
   deallocate(W_disk)
   !
   write (ch,'(a,i3.3)') "X_Q_",iq
   if (ver_is_gt_or_eq(ID,revision=308)) ch="X_Q_"//trim(intc(iq))
   !
   call io_bulk(ID,VAR=trim(ch),VAR_SZ=(/X%ng_db,X%ng_db,2,Xw%n(1)/))
   !
   do i1=1,Xw%n(1)
     !
     if (write_is_on(ID)) call mat_c2r(X_mat(:,:,i1),X_disk)
     !
     call io_bulk(ID,R3=X_disk,IPOS=(/1,1,1,i1/))
     !
     if (read_is_on(ID)) then
       if (different_db_RL_order)  then
         call mat_r2c_ordered(X_disk,X_mat(:,:,i1),X_RLcomp_ig)
       else
         call mat_r2c(X_disk,X_mat(:,:,i1))
       endif
     endif
     !
   enddo
   !
   deallocate(X_disk)
 endif
 !
 call Fragments_Restart(ID,current_fragment=iq,fragments_todo=nqibz)
 !
 ! No I/O errors I update ioX at the first q TODO
 !
 ioX=io_restart_point(ID)
 !
1 call io_disconnect(ID)
 !
 contains
   !
   integer function X_RLcomp_table()
   !--------------------------------
   !
   ! X has been written using a different RL vectors ordering.
   ! X_RLcomp_ig is defined in such a way that
   !
   ! X_( X_RLcomp_ig(i), X_RLcomp_ig(j) ) = X_on_disk_(i,j)
   !
   use vec_operate,   ONLY:v_is_zero
   use R_lattice,     ONLY:ng_in_shell,n_g_shells
   use zeros,         ONLY:define_zeros
   !
   integer :: is,ig1,ig2,g_limits(2)
   real(SP):: v(3),local_zero(3)
   !
   X_RLcomp_table=0
   X_RLcomp_ig=0
   !
   do is=1,n_g_shells
    if (ng_in_shell(is) > X%ng) then
      if (any(X_RLcomp_ig==0)) X_RLcomp_table=-1
      return
    endif
    !
    g_limits=1
    if (is > 1 ) g_limits = (/ ng_in_shell(is-1), ng_in_shell(is) /)
    !
    call define_zeros(vector_=g_vec,zero_=local_zero)
    !
    g_loop: do ig1=g_limits(1),g_limits(2)
      do ig2=g_limits(1),g_limits(2)
        !
        v=RL_vecs_disk(ig1,:)-g_vec(ig2,:)
        !
        if (v_is_zero(v,zero_=local_zero)) then
          X_RLcomp_ig(ig1)=ig2
          cycle g_loop
        endif
        !
      enddo
    enddo g_loop
   enddo
   !
   end function
   !
end function
